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INTRODUCTION 






T'tnis»study®liiassbee|i^focused on-Jhe development, validation and 
application of a fractional step solution method of the time- 
dependent incompressible Navier-Stokes equations in generalized 
coordinate systems/ A solution method that combines a finite- 
volume discretization with a novel choice of the dependent variables 
and a fractional step splitting to obtain accurate solutions in 
arbitrary geometries has-been previously developed for fixed-grids^ 

see-Ref r- ftv* hJ 


In the present research effort, this solution method is 
extended to include more general situations, including cases with 
moving grids. The numerical techniques are enhanced to gain 
efficiency and generality. This^report summarizes briefly the work 
performed during the period October 1, 1988 through February 15, 
1990. Additional details on the various aspects of the study are 
given in Appendix A.^ 

NUMERICAL ENHANCEMENTS 


The fixed grid solution method has been extended to general 
moving grids. Errors originating from the discrete approximation of 
the time-dependent coordinate system are minimized by satisfying 
the discrete geometric conservation laws for the time-varying 
computational cells. To improve the efficiency of the fixed grid 
case, two versions of the solution method have been coded: (1) 

fixed-grid method, (2) moving-grid method. 

During the present study, several enhancements of the 
numerical method have been introduced. A partial list of the 

modifications is given below: 

(1) Implementation of more general boundary conditions 
implicitly. The allowable boundary conditions are: 
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(a) Periodic conditions, 

(b) Symmetric conditions, 

(c) Mixed Dirichlet and Neumann type boundary conditions. 

(2) The solution method has been extended to geometrically 
singular boundaries for the three types of grid topologies (C, O 
or H grids). 

(3) A multi-grid Poisson solver has been written and partially 
debugged. 

(4) Extensive efforts have been made to increase the efficiency of 
the method by improved vectorization. Presently, the fixed- 
grid method runs at 80 MFLOPS on the CRAY YMP (single CPU) 
and about 300 - 400 . 10’ 6 CPU sec/mesh-point/time-step are 
consumed. The moving-grid code is not yet fully vectorized. 

VALIDATION OF THE METHOD 

Several additional cases have been solved to validate the 
method against other numerical and experimental results. In all the 
cases tested so far, good agreement is obtained. The validation 
cases include: 


Fixed-Grid Case: 

(1) Flow in a two-dimensional polar cavity. 

(2) Flow in a two-dimensional channel with a fixed constriction 
and a time variable pressure gradient. 

(3) Flow over a two-dimensional elliptic airfoil with a steady and 
pulsatile upstream flow and a high laminar Reynolds number 
(Re = 14,300). 

(4) Three-dimensional flow in curved ducts, both with rectangular 
and circular cross-sections. 

(5) Flow over a submarine body at an incidence of 0° and 20°. 
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Additional details of these validation cases can be found in 
Ref. (2). A brief summary of preliminary results for two 
validation cases is given below: 

Flow Over an Elliptical Airfoil at a High Reynolds Number 

The two-dimensional flow over an elliptical airfoil of 
thickness ratio 1:2.91 at 14° angle of attack and a Reynolds number 
of 14,300 has been solved to compare with the recent experimental 
results. Two cases have been considered, (a) steady upstream flow, 
(b) pulsatile upstream flow. In the second case, a sinusoidally 
pusatile upstream flow with an amplitude of 5% of the steady part, 
and a non-dimensional period of T = 6.86 was simulated. A non- 
orthogonal 0-type grid of 161 x 141 mesh points in the radial and 
circumferential directions, respectively, has been used. Figure 1 
gives the time-evolution of the lift and drag coefficients for the 
steady and pulsatile upstream flows (it should be noted that in the 
pulsating case, the time is normalized by the period time T). The 
analysis of the results and the comparisons with the experimental 
results will be reported elsewhere. 

Flow Over a Submarine Body at Low Reynolds Numbers 

The axisymmetric flow over a submarine body has been 
computed for a low Reynolds number of Re = 1000 and a zero angle of 
attack. Figure 2 compares the pressure coefficient on the body of 
the submarine with the computed results of Ryan (Private 
Communication), while Fig. 3 shows the effect of the Reynolds 
number on the pressure coefficient. Figure 4 shows the distribution 
of the pressure coefficient for an incidence of 20°- Figure 5 plots, 
for the same case, the limiting streamlines (viewed from the rear 
end of the submarine) and the particle traces (side view). 

Moving-Grid Case: 

(1) Flow over a circular cylinder with a moving outer boundary. 

(2) Flow in a two-dimensional channel with a moving constriction. 

(3) Flow in a two-dimensional cavity with a moving piston. 
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Additional details on these validation cases can be found in 
Appendix A. 

CONCLUDING REMARKS 

In the present study, a fractional step solution method of the 
time-dependent, viscous and incompressible Navier-Stokes 
equations has been extended, enhanced and validated for both fixed 
and moving generalized coordinate systems. The method has been 
used to simulate time-periodic vortical flow fields. Investigation 
of these flows by novel analysis methods that are being developed by 
the author, may advance the understanding of the complex vortical 
flow phenomena found in pulsating flows. 

The study has demonstrated the capabilities of the present 
fractional solution method in simulating accurately complicated 
incompressible time-dependent viscous flow fields. 
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Fig. 2 - Distribution of the pressure-coefficient on the submarine 
body at 0° incidence. 



Fig. 3 - Effect of Reynolds number on the pressure-coefficient on 
submarine body at 0° incidence. 
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Fig. 4 - Distribution of the pressure-coefficient on the submarine 
body at 20° incidence. 
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Fig. 5 - Limiting streamlines and particle traces for the submarine 
body at 20° incidence. 
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Abstract 

A solution method of the time- accurate, incompressible 
Navier-Stokes equations in generalised curvilinear moving 
coordinate systems is presented in this paper. Accuracy 
is achieved by a conservative finite-volume discretisation 
which satisfies the geometric conservation laws in gener- 
alised moving coordinate systems. The solution method 
is second-order accurate in space and first-order accurate 
in time. A fractional step solution method is used to 
efficiently solve the discrete equations. The unknowns, 
namely the pressure and the volume- fluxes, are chosen to 
facilitate the formulation of a consistent Poisson equa- 
tion and to obtain a robust Poisson solver with favorable 
convergence properties. The method is validated by com- 
parisons to other numerical and experimental solutions. 
The comparisons show good agreement. 

1 Introduction 

Numerous solutions have already been obtained for 
steady, incompressible flows in complicated three- 
dimensional configurations, see for example Kwak et a/. 1 . 
Time-dependent, viscous simulations require large com- 
putational resources and only recently, with the advent of 
new and large supercomputers, has their solution become 
feasible for complicated flow problems. However, the ma- 
jority of the existing studies consider cases where a fixed 
grid can be employed. In many realistic situations one or 
more boundaries move and therefore the choice of a fixed 
grid is inadequate. Moving boundaries can be found in 
internal flows with a moving piston, bio-fluid flows with 
elastic boundaries as in cardiac flows, as well as in exter- 
nal aerodynamics where the contour of the configuration 
may deform either due to forced motion of the controls 
or aero-elastic interactions. Another circumstance which 
requires the use of moving grids is the application of adap- 
tive grids to the solution of time-dependent flowfields. 

'Senior Scientist, MCAT Institute, Member AIAA 

(Research Scientist, Member AIAA 
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Very few flow solutions exist for cases with moving 
grids. Ogawa and Ishiguro 3 solved a two-dimensional 
model of the blood flow in a human ventricle and the 
dynamic stall process on an oscillating airfoil by employ- 
ing moving grids in a stream-function- vortidty formula- 
tion. Tamura et al * have solved the two-dimensional 
flow over a circular cylinder with forced vibrations and 
with vortex-induced vibrations using primitive variables. 
Rogers and Kwak 4 have recently extended their time- 
accurate high-order upwind solution procedure, based on 
the artificial compressibility method (see Ref. [5]), to 
moving grids. This method was applied to the solu- 
tion of the three-dimensional flowfield in a model of the 
artificial-heart developed by Pennsylvania State Univer- 
sity. Kiris and Rogers 6 used this method to solve the 
internal flow in a cavity with a periodically moving pis- 
ton. A somewhat restricted application of moving grids 
was used by Ralph and Pedley 7 to solve the flow in a 
two-dimensional channel with a moving indentation. The 
numerical method uses the stream-function-vorticity for- 
mulation and employs a specially devised time-dependent 
coordinate transformation to resolve the difficulties asso- 
ciated with the moving boundary. However, this proce- 
dure can be used only for a particular class of problems 
and can not be applied for more general cases where the 
mesh points may move arbitrarily. 

In all these works a finite-difference approach is used. 
The grid’s motion is accounted for by terms derived from 
the time-derivatives of the coordinate system. The finite 
volume approach to the discretisation of the incompress- 
ible Navier-Stokes equations offers flexibility in applying 
certain geometric conservation laws and minimising the 
truncation errors in complicated geometries. In the finite 
volume discretization methodology the geometric param- 
eters have clear meaning such as the volume and faces’ 
area of the computational cells. 

In the present work, a fractional step solution proce- 
dure is developed for the time-dependent, incompressible, 
viscous flows in generalized moving coordinate systems. 
The formulation of the governing equations, the dis- 
cretization process and the numerical solution phases are 
combined together to yield an accurate solution method 



for complies ted flow problems. Special attention is given 
to the satisfaction of the. "geometric conservation laws" 
(both in space and time) to minimise the discretisation 
errors in complicated coordinate systems. This work em- 
phasises the development of a conservative finite-volume 
formulation for moving grids with an efficient fractional 
step solution method based on a robust Poisson solver. 
Details pertaining to the fixed grid solution method can 
be found in other publications by the present authors. 8,0 

2 Formulation 

The equations governing the flow of isothermal, constant 
density incompressible, viscous fluids in a time-dependent 
control volume with the face S(i) and volume V(t) axe the 
conservation of mass 

^ dS • (u - v) = 0 , (1) 

and the conservation of momentum 

^jTudr = jf dS-t , (2) 

where t is the time, u is the velocity vector, dS is a surface 
area element and dV is a volume element. The surface 
element velocity resulting from the motion of the grid is 
v. The tensor T is given by 

t = -(u - v)u - P I + 1 / (Vu + ( Vu) T ) (3) 


3 Discretization 

3.1 Geometric Quantities 

A general nonorthogonal coordinate system ((, q,C) is de- 
fined (discretely) by 

* = r({,tj,C.<) * ( 6 ) 

where r = (x,y ,x) T is the Cartesian coordinate system 
and t is the time. The computational domain (£,q, C) i» 
divided into uniform primary cells with mesh sise = 
At) = AC = 1, and the center of each primary cell is 
designated by the indices i, j, n. The area of the face l of 
a primary cell, see Fig. ??, is given by the vector quantity 

dr dr . . 

s ~ d(7Ti) x d(T+2) ’ 

where the computational coordinates l = £, q or ( are in 
cyclic order and x is the cross product operator. The vec- 
tor S' has the magnitude of the face area and a direction 
normal to it. It is also the contravariant base vector VI 
scaled by the inverse of the Jacobian 1/7, i.e., S* — —VI, 

10 
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Vinokur 10 has pointed out that an accurate discretisa- 
tion should satisfy certain geometric identities. The con- 
dition that a cell is dosed should be satisfied exactly in 
the discrete form 

E s ' = 0 > (») 

i 


for Newtonian fluids. The quantity / is the identity ten- 
sor, Vu is the gradient of u while (-) T is the transpose 
operator. The pressure is P, and v is the effective kine- 
matic viscosity. 

The only differences between the fixed and the mov- 
ing grid equations are the terms that indude the surface 
dement vdodty v and the time dependence of the cell 
geometry (volume and face). The volume conservation of 
each time-varying cell (special case of (1)) requires 

^-jfdS-v = 0. (4) 

where the term dS • v represents the volume swept out by 
the face dS over the time increment dt. Thus, the mass 
conservation equation has exactly the same form as for 
fixed grids 

£ dS • u = 0 , (5) 

The usual practise is to transform equations (2) and (5) 
into a differential form. In the present work the integral 
formulation is maintained to assist in the derivation of 
the finite-volume equations for obtaining a conservative 
scheme for arbitrarily moving geometries. 


where the summation (with proper signs) is over all the 
faces of the computational cell. Equation (8) can be sat- 
isfied if S' is approximated from (7) by a proper approxi- 

&T 

mation of — (see details in Ref. [8]). The volumes of all 
discrete computational cells will sum up to the total vol- 
ume at a given time if the volume of each computational 
cell is computed by dividing the cell into three pyramids 
having in common the main diagonal and one vertex of 
the cell, resulting in the following 

v = 5 ( s ‘-i + s *-t + s *-i) 

' ( r »+}7+J.«+i ” 

In the present method, the volume conservation equa- 
tion (4) is satisfied discretdy by interpreting the term 
dS • v in equation (4) as the rate of the volume swept by 
the face dS. For example the volume swept by the face 
can be computed by a formula similar to (9) 

= j ((s<_,)‘ +«;_,+<_,) • 
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where the time-level is gives by k. The quantities 
and j are the areas swept by the motion of the face 
S< <-i - see the shaded areas in Fig. ??. The area 
is computed from 


«;.) = 5 (rtti -<+»> * _ tJ _, 

(u) 

and j can be computed in a similar way. The volume 
of the cell at the time level k + 1 is computed from (4) 


r *+i _ V k + £($y»)*+i , (12) 

l 


m = s* • u , (H) 

U< = S< • u , 

as the unknowns instead of the Cartesian velocity com- 
ponents. The quantities U ( , XT' and U ( are the volume- 
fluxes over the (, tj and C faces of a primary cell, respec- 
tively. In tensor algebra nomenclature, these are the con- 
travariant components of the velocity vector (in a stag- 
gered grid) scaled by the inverse of the Jacobian (1//). 
With this choice of the dependent variables, the continu- 
ity equation takes a form identical to the Cartesian case 

*&» - OJL, + u; ti - SJL, + I?*, - B<_ , = 0 . (IS) 


where the summation (with the proper signs) is over all 
the faces of the computational cell. Mote that 6V 1 / At 
is the volume-flux due to the motion of the coordinate 
system and has a meaning similar to the volume-flux U l . 
The computation of the volume of each cell at the time 
level k + 1 as the sum of the volume at time level k and 
the contribution of the volumes 6V 1 swept out by the cell's 
faces is important for the accurate representation of the 
momentum equation. 

In the present finite volume formulation, no coordinate 
derivatives appear directly in the discrete equations, as 
in the case of finite-difference formulas. Instead, quan- 
tities with clear geometric meaning such as the volume 
and faces area of the computational cells are used. Their 
discrete approximation is based on geometric interpreta- 
tions which satisfy the geometric conservation laws. A 
principal difference between the finite volume and the fi- 
nite difference approach to moving grids is in the inter- 
pretation of the quantity 6V 1 /At. In the finite- volume 
method it is treated as a geometric quantity which ex- 
presses the rate of displacement of a cell face, whereas in 
the finite-difference method, the grid velocity is combined 
with the fluid velocity to define a 'relative flow velocity' 
(see Ref. [10]). 

3.2 Mass Conservation Equation 


Accumulated experience with fractional step solution 
methods shows that the exact satisfaction of the discrete 
mass conservation equation is crucial for obtaining ac- 
curate solutions and for the convergence of the Poisson 
equation, see for example Ref [11]. Therefore, the simple 
form of (15), which can be satisfied to round-off errors in 
any coordinate system, suggests that the volume-fluxes 
are the ‘natural’ dependent variables for fractional step 
methods. This choice complicates the discretisation of 
the momentum equations, but is important for obtaining 
a divergence-free velocity field in generalised coordinate 
systems. 

3.3 Momentum conservation equation 

Spatial discretisation of the momentum conservation law 
(2) for a computational cell with volume V yields 

^(Vu) = ^ sI .J = F. 

Following the development of the scheme given by Rosen- 
feld et al. 8 , a general two-step temporal discretisation of 
(16) for a constant time-step At is 

(1 + e)(Vu)‘ +l - (1 + 2e)(Vu)‘ + efVu)*' 1 

= At(0F* +1 + (1 - 0)F* + *(F* - F*- 1 )), (17) 



The mass conservation equation (5) has the same form as 
for fixed grids. Second-order spatial discretisation over 
the faces of the primary computational cells yields 

(S* • n) j+ j - (S* • u )i_i + (S' • n) J+ j - (S* • u),_j 

+ (S<-u). +J -(S<.«)._ i =0 (IS) 

Note that throughout the present paper the default in- 
dices (»', j, n) and the default time-level (k + 1) are usu- 
ally omitted for simplicity. Each term on the left hand 
ride of (13) approximates the volume-flux over a face of 
the primary cell. A simple discretised mass conservation 
equation can be obtained by using the following variables 

U< = S«-n, 


where k is the time level and e, 6 and «j> are arbitrary pa- 
rameters. A truncation error analysis shows that second- 
order accuracy in time is obtained only if 

2 + * = $ + ^ • (18) 

In order to replace u by the new dependent variables U l , 
the corresponding area vectors are dotted with the mo- 
mentum equations. The integral momentum equation is 
applied on different computational cells for each unknown 
U l . Each cell has the dimensions x Atj x AC, but 

the centers are located at (* 4- (i,j 4- 1,») and 

(»,y',u + l) for the U*,U*, and U* momentum equations, 
respectively. The derivation of the (-momentum equation 
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will be described in this section. The other two momen- 
tum equations can be obtained by cyclic permutation. 

The dependent variables are the pressure and the 
volume-fluxes (see 3.2). In order to replace the unknown 
u by U l , eq. (17) is dotted by S* and using the identity 

\i = s ( u ( + s v ir i + s < v< = s m tr* (is) 

results in 

(1 + e)(V J7*)* +l - (1 + 2t)(V U') k ( St ■ S,‘) 

+ c(VP , )‘- , (S‘-S}- 1 ) (20) 

= At(0S< • F‘ +1 + (1 - fl)S f - F* + • (** - F*" 1 )) 

where S m is the inverse base of S 1 . Note that the k + 1 
is the default time-level and therefore S* is computed 
at k + 1. In order to save computations, only the term 
L = S* • F at the time-levels k — 1, k, and k + 1 should 
appear in (20). Therefore, the right hand side of (20) is 
modified to 

flg5( 20 ) _ ^ s< F} * +1 + (1 _ g ^ g( 

+ ^(S<-F)‘-(S<-F) k - 1 ) (21) 

= 0L k+1 + (1 - 0)L k + 4(L k - L k ~ l ) . 

This approximation is first-order in time. The tem- 
poral truncation errors can be minimised by choosing 
c = | and 4>— 1 — 0- The resulting scheme is 

i(V U<) k+1 - 4(V py ( s ( • sf ) + ( V U') k ~ l ( s ( • s? _1 ) 

= 2A t(0L k+l + (1 - 6)(2L k - L k ~ 1 )) . (22) 

For stability reasons, usually 0 = 1 is chosen. 

The computation of the operator L is similar to the 
fixed grid case and details can be found in 8,0 . The only 
difference is in the computation of the convection terms 
which should include the motion of the grid. For example, 
the convection flux of the {-momentum equation on the 
{-face center (*, j, n) is given by 



The difference equations could have been second-order 
accurate in time if dV* would not lag in time by — - over 

the volume-flux terms U l (see eq. (10)). 

The resulting discrete equations are conservative in any 
moving coordinate system and are spatially second-order 
accurate. For high Reynolds number flows, fourth-order 
dissipation is added to annihilate high frequency compo- 
nents of the solution. The dissipation terms are inter- 
preted in terms of fluxes and therefore the conservation 
properties of the equations are mmntMawi. 


4 Solution Method 

The dependent variables are the pressure, defined at the 
center of the primary cells, and the volume-fluxes de- 
fined on the faces of the primary cells. This selection is 
equivalent to a finite-difference formulation over a stag- 
gered grid with the choice of scaled contravariant velocity 
components as the unknowns. The mass and momentum 
conservation equations (5) and (22) are solved by a frac- 
tional step method. First the momentum equations are 
solved for an approximate U l by an approximate factor- 
isation method. In the second stage, the pressure and 
the volume-fluxes U l are corrected to satisfy the mass 
conservation equation by solving a Poisson equation. It 
should be noted that the exact satisfaction of the dis- 
crete mass conservation equation ensures the convergence 
of the Poisson equation with good convergence proper- 
ties. More details on the solution method can be found 
in Refs. [8,9]. 

The original solution method used an explicit approx- 
imation of the convection terms. In the present imple- 
mentation, the most important convection terms may be 
approximated also implicitly and therefore the scheme has 
no severe restrictions on the allowable CFL-number. Ex- 
plicit or implicit fourth-order dissipation can be used for 
high Reynolds number flows. 

5 Results 

Several test cases have been solved so far. The solution 
procedure and the computer-code are capable of solving 
three-dimensional fiowfields, but at this stage of the study 
only two-dimensional cases have been considered. 

5.1 Flow over a Circular Cylinder 

The fiowfield over a circular cylinder at Re = 40 was 
solved as the first test case. This example is a common 
test case for fixed grid solutions since the resulting flow- 
field is quite complicated and numerous other numerical 
and experimental results are available for comparison. In 
the present solution, the grid is expanding radially due to 
a moving circular outer boundary 

= ( 23 ) 

where Rm— is ike distance of the outer boundary from 
the cylinder center (normalised with the cylinder’s diam- 
eter). In the present case the dependence of the grid on 
time is pre-determined. . However, the same procedure 
may be used to solve a time-dependent problem by an 
adaptive-grid method, if a suitable adaptation criterion 
would have been defined. 

A cylindrical coordinate system is used with mesh 
points clustered near the cylinder and in the wake region. 


- 4 



Figure ?? gives the grid and the instantaneous stream* 
lines for t = 1 and t = 8. Between the times shown, 
the radial mesh si se is almost doubled in sise. Neverthe- 
less, the time-evolution of the separation bubble, which 
is a critical test parameter, agrees well with previous nu- 
merical and experimental results, see Fig. ??. The full 
line in Fig. ?? describes the solution obtained in the this 
study while the dashed line gives the solution obtained by 
the fixed-grid version of the present method. The trian- 
gular symbols and the dotted line give the experimental 
results of Coutanceau and Bouard 13 and the numerical 
solution of Collins and Dennis, 13 respectively. All the 
results practically coincide, except the numerical results 
of 13 for t > 7. As was explained by Rosenfeld ti al .,* this 
deviation results mainly from the "wall effect.” Note es- 
pecially the good agreement with the experimental results 
which have been obtained for an equivalent outer bound- 
ary of about 14 diameters. This is about the same as the 
outer boundary distance in the present computations at 
t as 8. 

5.2 Channel Flow with an Asymmetric 
Oscillating Indentation 

This test case considers the flowfield in a two-dimensional 
channel with an oscillating constriction. This geometry 
is a model for the large amplitude self-excited oscillations 
that arise when fluid flows through a collapsible tube such 
as a vein. In a series of flow visualisations, Stephanoff et 
al. 1 * and Pedley and Stephanoff 11 have found in the core 
of the flow a train of waves, downstream of the oscillating 
indentation, and a double row of eddies along the walls 
of the channel. 

In the experimental apparatus the walls of the chan- 
nel were rigid except for an indentation of length 10 (the 
distance between the fixed parts of the walls is the ref- 
erence unit length). The indentation is made of a thick 
rubber membrane and is driven by a piston with a sinu- 
soidal motion with time. At the begining of each cycle, 
the indentation is flushed with the channel’s wall. The 
channel starts at a distance of 120 units upstream of the 
oscillating constriction and is 250 units long. The flow- 
field upstream of the indentation is fully developed and 
the flow is essentially laminar. 

In the present work a simplified model of the test ap- 
paratus is employed. The upstream boundary is at a dis- 
tance of 5 units from the oscillating constriction and the 
downstream boundary is at a distance of 30 units from 
the upstream boundary. An algebraic grid is generated at 
each time step with 31 x 251 x 3 points along the heigth, 
length and width of the channel, respectively. Points are 
clustered near the two walls and a non-uniform distribu- 
tion of points is used in the axial direction with cluster- 
ing in the region downstream of the indentation where 
the flow structure was found experimentally to be moat 


interesting. The mesh points are fixed in time except at 
the indentation region where mesh points stretch linearly 
with the the distance between the moving indentation 
and the opposite fixed wall. The maximum indentation is 
0.38 units. The grid for this position is shown in Fig. ??, 
where only every other point is plotted in the normal and 
axial directions and the vertical scale is twice as large as 
the axial scale. The shape of the indentation is approx- 
imated by a hyperbolic tangent, as has been suggested 
by Pedley and Stephanoff 13 and Ralph and Pedley. 7 The 
apparent asymmetry of the indentation is a result of the 
non-uniform grid used in the present case. It was found in 
Ref. [7], as well as in the computations performed in this 
study, that the flowfield downstream of the indentation is 
not sensitive to the upstream conditions. 

Fully developed flow is given as the initial condition and 
as the upstream boundary condition. The downstream 
condition imposes parallel flow and a uniform pressure 
gradient which is computed from the integral mass con- 
servation law. At the upper and lower walls, the velocity 
of the wall is specified (sero everywhere, except at the 
moving parts of the boundary). 

Figure ?? shows a comparison of the instantaneous 
streamlines with a flow visualisation of Stephanoff et al. 1 * 
at the non-dimensional time of t = 0.55 (based on the 
period) and for a case with a Strouhal number of St = 
0.038 and a Reynolds number of Re = 610. The flow 
visualisation shows quite a complicated flowfield although 
the resolution is not good enough for revealing the fine 
details, especially in the regions of eddy motion. The first 
separation of the flow occurs downstream of the sloping 
wall. A second large eddy is formed on the opposite wall 
with a secondary separation bubble buried inside it. Still 
further downstream, an additional pair of weaker vortices 
appear, one on each wall. The core of the flow prescribes 
a wave motion which can be observed easily both in the 
experimental and the numerical results. 

Favorable agreement can be obtained if the plot of the 
numerical result is moved about 0.4 units in the down- 
stream direction relative to the experimental results. In 
other words, the separation length of the first eddy at 
the upper-wall is under-predicted in the present compu- 
tations. Nevertheless, the distances between the vortices, 
which are related to the wavelength of the core-flow, com- 
pare favorably. The reason for this discrepancy is not 
clear and requires additional study. The differences do 
not seem to be related to the resolution since grid re- 
fining does not affect the separation length significantly. 
Kiris and Rogers 6 has solved the same case with a high- 
order upwind scheme and found similar results. The fact 
that two solution procedures, which differ in numerous 
aspects, result in rimflar solutions can indicate that the 
differences may be related to the non-exact reproduction 
of the experimental set-up and conditions. 

One very likely reason for the discrepancy may be at- 
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tributed to three-dimensional effects. In the present com- 
putations a two-dimensional flow is assumed. In the ex- 
periments conducted by Stephanoff e< al. 14 the width of 
the channel was 10 times larger than its height. The 
three-dimensional affects of the sidewalls on the center- 
plane, where the visualisations were taken, eonld not be 
assessed. Armaly et at. 19 have investigated experimen- 
tally the steady flow past a backward-facing step and 
found three-dimensional effects for Re > 400. The most 
sensitive quantity was the length of the separated eddies. 
In unsteady flows, the three-dimensional effects are even 
more pronounced, see for example Ref. [17]. Another rea- 
son for the disagreement could be the inaccurate descrip- 
tion of the wall shape in the numerical computations. Ob- 
viously, any difference in the shape may affect the length 
of the separation from the indentation. 

The flow evolution for Re = 600 and a higher Strouhal 
number, St = 0.057, is given in Fig. ??. Here the 
instantaneous streamlines are plotted for several non- 
dimensional times at the first cycle of the indentation’s 
motion. The flow development is essentially the same 
as that found experimentally by Pedley and Stephanoff 18 
and is summarised only briefly here. 

Near the beginning of the cycle, a single separation 
bubble forms on the sloping wall of the indentation 
(Fig. ??a) and the core-flow remains parallel to the walls. 
As time increases, the separation length increases and 
a second counter-rotating eddy appears on the opposite 
wall, downstream of the first eddy - see Fig. ??b. Still 
later, a third eddy is formed at the upper wall further 
downstream (Fig. ??c). As time proceeds, the core-flow 
becomes wavy and a sequence of eddies with alternating 
signs appear on the walls. The amplitude of the core- 
flow increases with time until t = 0.75 (Fig. ??g). In 
the last quarter of the period, the eddies shrink in sise 
and strength and are washed downstream. At the end of 
the cycle the residual eddies are quite small and do not 
affect the next cycle. We shall adopt the convention of 
Ref. [15] in labeling the eddies alphabetically, as shown 
in Fig. ??d. 

Of particular interest is the phenomenon of eddy- 
doubling which was found experimentally (Ref. [15]) as 
well as computationally (Ref. [7]). In this phenomenon, a 
single eddy splits into two co-rotating eddies. The present 
calculations show that eddy-doubling occurs for eddies A, 
B and C for t > 0.55. 

The genesis of the wavy flow as well as the mechanism 
of the eddy-doubling is still not fully understood. Most 
previous studies agree that the wavy core-flow is deter- 
mined primarily by in viscid vortidty dynamics and is a re- 
sult of the non-sero vortidty gradient, see Refs. [7,15,17]. 
It is also dear that the eddy-doubling is essentially a 
viscous bifurcation, Refs.[7,15]. The motivation of the 
present paper is to describe the numerical solution proce- 
dure, therefore no further attempt will be made to analyse 


this flowfidd. 

Figure ?? plots the center of the vortices A, B, C 
and D (which approximatdy coinddes with the crests 
and troughs of the core flow streamlines) as a function of 
time for four different results. The experimental results 
are given by square symbols while the various numeri- 
cal computations are represented by lines. The numer- 
ical results indude the present study (full lines)as well 
as the computations of Kiris and Rogers 4 (dash lines) 
and Ralph and Pedley 7 (dotted lines). The agreement 
between the present results and the results obtained by 
Kiris and Rogers is good, except perhaps for the first 
separated eddy (A). The numerical results of Ralph and 
Pedley, 7 who used a stream-function-vortidty formula- 
tion and a mesh with about 3.5 times more points in the 
axial direction, show better agreement with the experi- 
mental results, obviously due to the better prediction of 
the first separation length. However, the distances be- 
tween the vortices (or the wavdength) compare favorably 
between all these results. 

5.3 Internal Flow Driven by a Piston 

The last case demonstrates the capabilities of the present 
numerical procedure. In this case, a two-dimensional in- 
ternal flow in a rectangular cavity driven by a vertically 
pulsating piston is solved for a Reynolds number Re=100, 
see geometry in Fig. ??. During the downward stroke of 
the piston, the fluid leaves through an exit at the left ver- 
tical wall, and during the upward motion the entrance at 
the right vertical wall opens and the exit is dosed. This 
flowfidd is an idealisation of a reciprocal engine or an ar- 
tificial heart. The piston has a sinusoidal motion between 
y — 0.6 and y = 0.3 with a period of 1.2 non-dimensional 
seconds. 

An algebraic grid of 49 x 45 x 3 points in the z, y and z 
directions, respectivdy has been employed. Mesh points 
are dustered near the two side walls and near the lower 
wall. Along the y-direction, the mesh points are stretched 
linearly with the motion of the piston. The z-direction 
distribution remains fixed. Figure ?? shows the time 
evolution of the foree-coefiidents in the y and z direc- 
tions on the lower stationary horizontal wall (opposing 
the piston). A fully devdoped periodic flow is seen to 
persist from the second cyde of the piston. The abrupt 
dosing of the exit (at t = 0.6, 1.8, 3., etc.) and the abrupt 
dosing of the entrance (at t — 1.2, 2.4, 3.6, etc.) can be 
dearly observed as a sharp change in the x -component of 
the foree-coeffident. 

Figure ?? gives the vdodty vector and the pressure 
contours for the fifth cyde. The vdodty vector is {dot- 
ted for every other mesh point in both z and y directions. 
At the be ginning of the cyde a strong vortex, which was 
created in the previous cyde near the entrance, is found 
along with a contra-rotating weaker vortex. At half stroke 
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downward, no vortices exist (except very weak vortices 
near the two lower corners). A very strong pressure gradi- 
ent bnilds np at the exit. As the entrance opens, Fig. ??c, 
two main vortices are found in the core region, bat they 
are weaker than at the beginning of the cycle. Daring 
the upward motion of the piston, a strong vortex is gen- 
erated at the entrance, Fig. ??d. A favorable agreement 
has been foond with Kiris and Rogers, 6 who has solved 
the same case by a time-accurate artificial compressibility 
method with high-order upwinding. 

6 Concluding Remarks 

A general solution method of the time-accurate incom- 
pressible Navier-Stokes equations in a generalised curvi- 
linear moving coordinate system is presented. Accuracy 
is achieved by employing conservative finite-volume dis- 
cretisation with special attention to the satisfaction of 
geometric conservation laws in generalised moving coor- 
dinate systems. A fractional step solution method is used 
to efficiently solve the discrete equations. The pressure 
and the volume-fluxes are chosen to be the independent 
variables. This choice simplifies the construction of a ro- 
bust Poisson solver. The present method is an alternative 
to other finite difference methods. 

Several test solutions show promising results. It is es- 
pecially encouraging to obtain smooth solutions near in- 
flow and outflow boundaries which are time dependent. 
Future work will include the solution of more realistic 
three-dimensional cases with moving boundaries as well 
as time-dependent problems with self-adaptive grids. 
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2: Definitions related to moving grids 
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Figure 5: The grid for the channel flow with oscillating indentation 
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Figure 6: Comparison of the instantaneous streamlines at t = 0.55 and Si = 0.038, Re = 810. (a) exp erim ental 
results (b) present results 
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Figure 7: Instantaneous streamlines, St 
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Figure 11: Velocity vector and pressure contours 











